June 19, 2017

Get: dataRetrieval for NWIS/WQP data

library(dataRetrieval)
o2data_raw <- readNWISuv(
  siteNumbers='02336300', parameterCd='00300', 
  startDate='2013-01-01', endDate='2017-06-18', tz='UTC')
o2data <- renameNWISColumns(o2data_raw)

Get: dataRetrieval for NWIS/WQP data

##    agency_cd  site_no            dateTime DO_Inst DO_Inst_cd tz_cd
## 1       USGS 02336300 2013-01-03 05:00:00     9.3          A   UTC
## 2       USGS 02336300 2013-01-03 05:15:00     9.4          A   UTC
## 3       USGS 02336300 2013-01-03 05:30:00     9.4          A   UTC
## 4       USGS 02336300 2013-01-03 05:45:00     9.4          A   UTC
## 5       USGS 02336300 2013-01-03 06:00:00     9.4          A   UTC
## 6       USGS 02336300 2013-01-03 06:15:00     9.5          A   UTC
## 7       USGS 02336300 2013-01-03 06:30:00     9.5          A   UTC
## 8       USGS 02336300 2013-01-03 06:45:00     9.5          A   UTC
## 9       USGS 02336300 2013-01-03 07:00:00     9.5          A   UTC
## 10      USGS 02336300 2013-01-03 07:15:00     9.4          A   UTC

Get: geoknife for spatial data

library(geoknife)
stencil <- simplegeom(data.frame(c(-84.40764,33.82031)))
daymet <- 'https://thredds.daac.ornl.gov/thredds/dodsC/daymet-v3-agg/na.ncml'
fabric <- webdata(url=daymet, variable='srad', times=c('2013-01-01','2017-01-01'))
job_out <- geoknife(stencil, fabric, REQUIRE_FULL_COVERAGE='false', wait=TRUE)
gdpdata <- result(job_out, with.units = TRUE)
names(gdpdata)[2] <- 'srad'

Get: geoknife for spatial data

##      DateTime  srad variable statistic units
## 1  2013-01-01 224.0     srad      MEAN  W/m2
## 2  2013-01-02 124.8     srad      MEAN  W/m2
## 3  2013-01-03 195.2     srad      MEAN  W/m2
## 4  2013-01-04 278.4     srad      MEAN  W/m2
## 5  2013-01-05 297.6     srad      MEAN  W/m2
## 6  2013-01-06 236.8     srad      MEAN  W/m2
## 7  2013-01-07 297.6     srad      MEAN  W/m2
## 8  2013-01-08 284.8     srad      MEAN  W/m2
## 9  2013-01-09 214.4     srad      MEAN  W/m2
## 10 2013-01-10 198.4     srad      MEAN  W/m2

Explore: dygraphs for interaction

library(dygraphs)
library(xts)
o2data_xts <- xts(o2data$DO_Inst, o2data$dateTime)
d <- dygraph(o2data_xts)

Explore: dygraphs for interaction

Explore: ggplot for complexity

library(ggplot2)
g <- ggplot(o2data, aes(x=dateTime, y=DO_Inst, color=DO_Inst_cd)) +
  geom_line() + theme_bw()

Explore: ggplot for complexity

Clean: remove flagged values

library(dplyr)
o2data_accepted <- filter(o2data, DO_Inst_cd=='A')

Clean: remove flagged values

##    agency_cd  site_no            dateTime DO_Inst DO_Inst_cd tz_cd
## 1       USGS 02336300 2013-01-03 05:00:00     9.3          A   UTC
## 2       USGS 02336300 2013-01-03 05:15:00     9.4          A   UTC
## 3       USGS 02336300 2013-01-03 05:30:00     9.4          A   UTC
## 4       USGS 02336300 2013-01-03 05:45:00     9.4          A   UTC
## 5       USGS 02336300 2013-01-03 06:00:00     9.4          A   UTC
## 6       USGS 02336300 2013-01-03 06:15:00     9.5          A   UTC
## 7       USGS 02336300 2013-01-03 06:30:00     9.5          A   UTC
## 8       USGS 02336300 2013-01-03 06:45:00     9.5          A   UTC
## 9       USGS 02336300 2013-01-03 07:00:00     9.5          A   UTC
## 10      USGS 02336300 2013-01-03 07:15:00     9.4          A   UTC

Analyze: low-O2 events

A common O2 criterion: daily mean >= 5 mg/L, daily minimum >= 4 mg/L

library(dplyr)
o2data_daily <- o2data_accepted %>%
  mutate(date=as.Date(dateTime)) %>%
  group_by(date) %>%
  summarize(DO_Mean=mean(DO_Inst),
            DO_Min=min(DO_Inst),
            DO_Max=max(DO_Inst)) %>%
  mutate(Violation = DO_Mean < 5.0 | DO_Min < 4.0)

Analyze: low-O2 events

g <- ggplot(o2data_daily, aes(x=date)) +
  geom_ribbon(aes(ymin=DO_Min, ymax=DO_Max), fill='blue', alpha=0.2) +
  geom_line(aes(y=DO_Mean), color='blue') +
  geom_point(data=filter(o2data_daily, Violation), aes(y=DO_Mean), color='red') +
  theme_bw()

Analyze: low-O2 events

Analyze: low-O2 events

o2_violations <- filter(o2data_daily, Violation)
## # A tibble: 17 x 5
##          date  DO_Mean DO_Min DO_Max Violation
##        <date>    <dbl>  <dbl>  <dbl>     <lgl>
##  1 2014-06-21 5.411458    3.9    6.5      TRUE
##  2 2014-06-22 4.927174    3.7    5.9      TRUE
##  3 2014-09-19 5.992708    3.6    8.0      TRUE
##  4 2015-05-26 6.348958    3.1    9.1      TRUE
##  5 2015-05-27 6.955914    3.6    8.4      TRUE
##  6 2015-07-02 5.075789    3.6    6.5      TRUE
##  7 2015-07-23 4.981250    4.2    6.1      TRUE
##  8 2015-07-30 4.954167    4.3    5.8      TRUE
##  9 2015-07-31 4.833333    4.0    5.8      TRUE
## 10 2015-08-02 4.813542    4.3    5.6      TRUE
## 11 2015-08-03 4.733333    4.0    5.6      TRUE
## 12 2015-08-05 5.023656    3.9    5.9      TRUE
## 13 2015-08-06 4.427368    3.5    5.3      TRUE
## 14 2015-08-07 4.925000    4.1    5.4      TRUE
## 15 2015-08-09 4.820000    4.3    5.5      TRUE
## 16 2015-08-10 5.038947    3.8    6.5      TRUE
## 17 2016-10-22 5.426042    3.9    8.1      TRUE

Analyze: streamMetabolizer

nwisdata <- renameNWISColumns(readNWISuv(
  siteNumbers='02336300',
  parameterCd=c('00300','00060','00010'), # oxygen (mg/L), discharge (cfs), water temperature (C)
  startDate='2015-01-01', endDate='2016-01-01', tz='UTC'))
coords <- readNWISsite('02336300')

library(streamMetabolizer)
metabdata <- nwisdata %>%
  mutate(
    solar.time = convert_UTC_to_solartime(dateTime, longitude=coords$dec_long_va, time.type='mean solar'),
    discharge = Flow_Inst * 0.0283168, # cfs to cms
    depth = calc_depth(discharge), # uses Raymond et al. 2012 to estimate depth from discharge
    DO.sat = calc_DO_sat(Wtemp_Inst, calc_air_pressure(elevation=coords$alt_va * 0.3048)), # serious modeling would use better inputs here
    light = calc_light(solar.time, latitude=coords$dec_lat_va, longitude=coords$dec_long_va)) %>%
  select(solar.time, DO.obs=DO_Inst, DO.sat, depth, temp.water=Wtemp_Inst, light)
# mm <- metab(specs(mm_name('bayes')), data=metabdata)

Analyze: metabolism at scale (Powell Center synthesis)

library(mda.streams)
mm_bayes <- 
  get_ts(
    c('sitedate_calcLon','gpp_estBest','er_estBest','K600_estBest'),
    site_name='nwis_02336300') %>%
  unitted::v() %>%
  filter(sitedate >= '2013-01-01', sitedate <= '2017-01-01')

Analyze: metabolism at scale (Powell Center synthesis)

library(cowplot)
g <- plot_grid(
  ggplot(mm_bayes, aes(x=sitedate, y=gpp)) +
    geom_hline(yintercept=0, color='darkgrey') +
    geom_point(color='#007929', alpha=0.2) +
    theme_bw() + xlab('') + ylab(expression(GPP~(gO[2]~m^{-2}~d^{-1}))),
  ggplot(mm_bayes, aes(x=sitedate, y=er)) +
    geom_hline(yintercept=0, color='darkgrey') +
    geom_point(color='#A64B00', alpha=0.2) +
    theme_bw() + xlab('') + ylab(expression(ER~(gO[2]~m^{-2}~d^{-1}))),
  nrow=2, align='v')
## Warning: Removed 245 rows containing missing values (geom_point).

## Warning: Removed 245 rows containing missing values (geom_point).

Analyze: metabolism at scale (Powell Center synthesis)

Analyze: rolling means

mm_bayes_rollmean <- mm_bayes %>%
  filter(., complete.cases(.)) %>%
  mutate(
    gpp_roll = zoo::rollmean(gpp, k=15, na.pad=TRUE),
    er_roll = zoo::rollmean(er, k=15, na.pad=TRUE),
    K600_roll = zoo::rollmean(K600, k=15, na.pad=TRUE))

Analyze: rolling means

g <- plot_grid(
  ggplot(mm_bayes_rollmean, aes(x=sitedate, y=gpp)) +
    geom_hline(yintercept=0, color='darkgrey') +
    geom_point(color='#007929', alpha=0.2) +
    geom_line(aes(y=gpp_roll), color='#007929', na.rm=TRUE, size=1) +
    theme_bw() + xlab('') + ylab(expression(GPP~(gO[2]~m^{-2}~d^{-1}))),
  ggplot(mm_bayes_rollmean, aes(x=sitedate, y=er)) +
    geom_hline(yintercept=0, color='darkgrey') +
    geom_point(color='#A64B00', alpha=0.2) +
    geom_line(aes(y=er_roll), color='#A64B00', na.rm=TRUE, size=1) +
    theme_bw() + xlab('') + ylab(expression(ER~(gO[2]~m^{-2}~d^{-1}))),
  nrow=2, align='v')

Analyze: rolling means

Analyze: metabolism vs light

ggplot(gdpdata, aes(x=DateTime, y=srad)) + geom_point()

Analyze: metabolism vs light

lightmetab <- inner_join(mutate(gdpdata, date=as.Date(DateTime)), mm_bayes_rollmean, by=c(date='sitedate'))

Analyze: metabolism vs light

##    DateTime.x  srad variable statistic units       date
## 1  2013-01-04 278.4     srad      MEAN  W/m2 2013-01-04
## 2  2013-01-05 297.6     srad      MEAN  W/m2 2013-01-05
## 3  2013-01-07 297.6     srad      MEAN  W/m2 2013-01-07
## 4  2013-01-08 284.8     srad      MEAN  W/m2 2013-01-08
## 5  2013-01-09 214.4     srad      MEAN  W/m2 2013-01-09
## 6  2013-01-11 121.6     srad      MEAN  W/m2 2013-01-11
## 7  2013-01-12 153.6     srad      MEAN  W/m2 2013-01-12
## 8  2013-01-13 214.4     srad      MEAN  W/m2 2013-01-13
## 9  2013-01-18 240.0     srad      MEAN  W/m2 2013-01-18
## 10 2013-01-19 336.0     srad      MEAN  W/m2 2013-01-19
##             DateTime.y         gpp        er     K600  gpp_roll   er_roll
## 1  2013-01-04 17:36:42 -0.06001882 -5.350751 8.031714        NA        NA
## 2  2013-01-05 17:36:42  0.10317484 -6.242426 8.494305        NA        NA
## 3  2013-01-07 17:36:42  0.19793295 -4.061600 7.016444        NA        NA
## 4  2013-01-08 17:36:42  0.09559557 -4.012591 6.702252        NA        NA
## 5  2013-01-09 17:36:42  0.14706421 -4.081570 7.311908        NA        NA
## 6  2013-01-11 17:36:42  0.35984947 -4.669400 6.816117        NA        NA
## 7  2013-01-12 17:36:42  0.38414785 -5.166045 8.194363        NA        NA
## 8  2013-01-13 17:36:42  0.47046427 -3.718497 6.200358 0.2430295 -3.962039
## 9  2013-01-18 17:36:42  0.12745140 -3.583609 5.265111 0.2773623 -3.763658
## 10 2013-01-19 17:36:42  0.29362743 -4.549904 8.269185 0.2752185 -3.612694
##    K600_roll
## 1         NA
## 2         NA
## 3         NA
## 4         NA
## 5         NA
## 6         NA
## 7         NA
## 8   7.339964
## 9   7.297291
## 10  7.161015

Analyze: metabolism vs light

g <- ggplot(lightmetab, aes(x=srad, y=gpp)) +
  geom_point(na.rm=TRUE) +
  xlab(expression('Shortwave radiation'~(W~m^{-2}))) +
  ylab(expression('GPP'~(gO[2]~m^{-2}~d^{-1})))